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ABSTRACT 

A theoretical method for representing the physical characteristics 
of a portion of a mechanical piping system in terms of a stiffness 
matrix involving only those joints at which the sub-system is connected 
into the overall system is presented herein. The use of this method in 
reducing the size of the matrices that are needed for finding the normal 
vibration frequencies of a piping system is discussed, and an outline for 
a digital computer program to use this method is presented. 
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CHAPTER I 



INTRODUCTION 



1.1 General Remarks 

There are many aspects to be considered in the design of any pip- 
ing system; stresses, strains, pressure drops, and corrosion, to name 
a few. This thesis deals with only one aspect of the overall problem, 
that of the dynamic behaviour of the system as represented by its nat- 
ural frequencies of vibration, and presents a method whereby a large 
system may be systematically attacked by considering small sub-systems 
in order to improve computer utilization. 

Except for very simple or trivial problems the number of calcula- 
tions involved and the size of the matrices used make the use of a 
large digital computer mandatory. No digital computer program has been 
written, however, but an outlined plan of approach is included in this 
thesis. 

1.2 Scope of Work Presented 

Baird [1]* has presented a method for determining analytically the 
undamped natural frequencies and mode shapes of large three dimensional 
piping systems by use of stiffness matrices. He also has made use of 
topological matrices to express mathematically certain physical or 
boundary conditions. The most obvious inherent disadvantage of this 
method, however, is the large size of the matrices which may be re- 
quired with the result that supplementary memory of some sort must be 

★Numbers appearing in this manner refer to the bibliography, page 41. 
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used in the computer solutions, thus increasing the computation time 
considerably. 

A technique is presented in this thesis whereby small portions of 
the system, called ganglia or sub-systems, can be analyzed and reduced 
to some sort of pseudo-elements having certain calculated characteris- 
tics with the result that only the reactions at the ends of these pseudo 
elements are of interest, and all the internal characteristics are en- 
gulfed in the general stiffness matrix for each ganglion. These individ 
ual sub-systems can be used to build up an entire system whose vibra- 
tional characteristics can be represented far more compactly than can 
those of the unsimplified system. 

An outline for a proposed computer program to utilize this method 
is included in Chapter IV. 

1.3 Definitions 

A number of definitions must be established and the reader is warn- 
ed that they are somewhat at variance with those of Baird [1J. The 
nomenclature used (see section 1.4) will also be at variance with Baird' 
Element - A piece of the system having two, and only two, ends and 
no internal branches. Normally an element will be a simple length 
of pipe. A tee, though, which is a more complicated item, should 
be thought of as three elements joined at a common point. Other 
items of piping hardware may be analyzed similarly. 

Universe - This term is used in the same sense in which it is used 
in statistics, namely the larger body of effects from which the 
system is taken. In this case the universe for a piping system 
consists of all the boundary forces and constraints, foundations, 
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etc. which affect the system. When a sub-system is isolated, its 
universe is the larger system itself. 

Node - A terminal point of an element is called a node. It is 
noted that a node is a theoretical point and not a physical item. 
There must be at least one element end at each node, though there 
may be any number greater than one. Nodes may be further classi- 
fied as: 

1. Trivial nodes - Those nodes upon which only two element 
ends are incident. 

2. Boundary Nodes - Those nodes at which there are forces and 
constraints applied to the system by the universe. For pur- 
poses of convenience when dealing with coordinate transforma- 
tions it is required that only one element be incident upon a 
boundary node. In many cases this may require the introduc- 
tion of an element connecting a primary node to the boundary 
node at which the system connects to the universe. These 
primary and boundary nodes may have the same geometrical loca- 
tion, though they are topologically distinct. 

3. Primary Nodes - All nodes of the system that are not trivial 
or boundary nodes. 

Path - A sequence of elements connecting two nodes such that there 
are only trivial nodes connecting these elements. 

Primary Path - A path connecting two primary nodes or a primary node 
and a boundary node. A primary path may contain only trivial nodes 
on its interior. 

Primary Element - An element made up of all the elements on a primary 



path. 
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State Vector - A column matrix of displacements and corresponding 



forces at any point. Since in general it requires three translations 
and three rotations to describe motion completely the state vector 
is a 12 x 1 matrix. For the purposes of this thesis it has been 
convenient to make the first six elements of the column the displace- 
ments in the six independent degrees of freedom, and the last six 
elements the corresponding forces. 

Transfer Matrix - A square (12 x 12) matrix representing the physi- 
cal characteristics of an element such that the state vector at one 
end will equal the product of the transfer matrix and the state 
vector at the other end, in that order. 

Dynamic Transfer Matrix - A transfer matrix that takes into account 
the mass of the elements and the frequency of vibration in addi- 
tion to the static properties (i.e., Young's modulus, Poisson's 
ratio, second moments of area, etc.). 

Ganglion - A sub-system "removed" from the system at primary nodes. 
"Removed" is used to indicate that the ganglion is considered as a 
separate system whose boundary nodes are the nodes at which it was 
connected into the overall system. 

1.4 Nomenclature 

A capital letter will normally be used to represent a matrix and 
appropriate sub- and super-scripts may be used to further define it as 
follows: 




a - number of rows 
b - number of columns 
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c - may be 1) blank, 2) T indicating the transpose of the 



matrix, or 3) -1 indicating the inverse of the matrix 



d - used for identification and described as required in 



context . 



Brackets, [ ], may be used to delineate a matrix or the partition of a 
matrix if the meaning is thereby clarified. They will always be used 
when the elements of a matrix are exhibited in extensio . Braces, J, 
will be used to permit the writing of the terms of a column vector as 



terms inside the braces are the n diagonal terms of an otherwise null 
matrix, though the terms themselves may represent square sub-matrices. 

In no case will more than one capital letter be used for a matrix; 
thus the expression AB will mean the product of matrices A and B in that 
order and will imply that A is conformable to B for multiplication. One 
should recall that matrix multiplication is distributive and associative 
though not necessarily commutative. 

All nodes will be numbered, and for convenience the boundary nodes 
will be given the lowest numbers. The ends of the elements will be 
identified by lower case English letters, two per element. 

Notation generally follows that of Pestel and Leckie [2]. The 
following are the matrix symbols used. 



A Node incidence matrix 

D, F Displacement and force vectors (or partitions) respectively 
G Coordinate rotation matrix 
K Dynamic stiffness matrix 
Node stiffness Matrix 



a row vector to save space. 




mean that the n 




Frequency determinant 
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K Primitive stiffness matrix 
P 

T Transformation matrix 
U Transfer matrix 
Z State vector 

By the use of subscripts any of the above symbols may be used to 

represent an element or partition of the full matrix, thus Z. - would 

1 , Z 

be an element or partition (depending upon the context) of matrix Z, 
namely the one appearing in row 1, column 2. 

The following scalar quantities are used: 

B The number of boundary nodes of a system 
M/2 The number of primary paths in a system 

(Therefore M is the number of primary element ends in the 
system.) 

N The number of nodes (boundary and primary) in a system 
The use of a prime, ('), will generally indicate that the quantity desig- 
nated refers to a sub -system. 
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CHAPTER II 



THE TRANSFER AND STIFFNESS MATRIX SOLUTION 



A brief summary of the solution by use of stiffness matrices for 
the characteristic (normal) frequencies of a piping system with general 
topology is given in this chapter in order to refresh the reader's memory. 
For further details of the theory see Baird [1], but again, be aware of 
notational differences. 

Before we discuss the dynamic stiffness matrix approach as develop- 
ed in [1] and [2] it is necessary to remind the reader that this direct 
mathematical approach deviates somewhat from the approach to planning a 
digital computer program because, although on paper we can manipulate 
symbols, in the computer we must deal with numbers. It is assumed that 
the user of this method has available the necessary dynamic transfer 
matrices for whatever fundamental elements are involved and can compute 
(in one way or another) the geometrical matrices necessary to relate all 
quantities to some single coordinate system. 



a. Calculate a transfer matrix for each primary element, eliminat- 
ing the trivial nodes along the path. Thus: 

U, - U CD.fi, ...D. 0 G.U, . 
k,m n,m n n-l,n n-1 1,2 1 k,l 

where k and m are element ends at primary and/or boundary nodes, and n. 



n-1, ...1 are the trivial nodes along the path. G , G ,,...G, are the 

n n-1 1 

coordinate rotation matrices providing the necessary coordinate align- 
ment at each node. The resultant transfer matrix, U, is in the coor- 

k,m 

dinate system of the element end at the node at which the path starts, 
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end k in this case. The relationship between k and m is now: 



Z = U. Z, 
m k,m k 



2.1 



b. Transform the transfer matrix of each primary element into a 
dynamic stiffness matrix as indicated in what follows; The user is cau- 
tioned to establish and adhere to an order within his state vector. 

The order used by Baird [1] is recommended. Referenced to his local 
coordinate system it is: 

Displacement in an x direction 
Displacement in a y direction 
Displacement in a z direction 
Rotation about an x axis 
Rotation about a y axis 
Rotation about a z axis 
Force in an x direction 
Force in a y direction 
Force in a z direction 
Moment about an x axis 
Moment about a y axis 
Moment about a z axis 

Partitioning the state vectors into force and displacement vectors, 
and the dynamic transfer matrix into four 6x6 matrices yields: 



m 



m 



as 


U 11 


u i7 




_ U 21 


U 22__ 



2.2 



An explanation of the use of the lower case f will follow shortly. The 
above equations then can be written: 
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m 



U ll D k + U 12 f k 



m 



U 21 D k + U 22 f k 



Solving the first of these for f^ gives 

£ k ’ "n D . ' »i'2 D n D k 



and substituting this into the second gives 



m 



D 21 D k + D 22< D l'2 D n - "ll V 



Reassembling these last two equations in matrix form, we have 

■1 



— — 




-f 


m 


k 




F 




m 




— . 





°n 



°2r 0 22 U 12 U ll 



■°12 

U 22 U l2 



m 



2.3 



There is a sign convention that must be carefully delineated at this 

time. In equation 2.1 the force elements of the state vector, Z^, (shown 

in 2.2 as the partitioned sub-matrix f^) are the forces applied to 

element km by the node at k, whereas the force elements of the state 

vector Z are the forces applied by element km to the node at m. It 
m — 

is obvious then that the negatives of the elements of sub-matrix f^ 
represent the forces applied b£ element km to the node at k; and if we 
define a new vector F^ = -f^* then both and F^ represent forces ap- 
plied by the element to the nodes. ^Then equation 2.3 can be rewritten: 



h 


B 

* 

II 


i 


F 




D 


L- raJ 




<— m_l 



where Ym> 

matrix in 2 



the dynamic stiffness matrix for the element km, 
3. 



is the square 
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c. The dynamic stiffness matrix of each primary element is now 
known in the coordinate system of the element end at which calculations 
started. In order to relate all these matrices to each other it is 
necessary that some global coordinate system and reference axes be 
chosen and that all the dynamic stiffness matrices be transformed into 
this system. Detailed discussion of the necessary geometric matrices 
and multiplications may be found in [1] and [3] and will not be repeated 
here. 

d. The system now consists of: 

M/2 primary nodes 
B boundary nodes 
N nodes (boundary and primary) 
zero trivial nodes 



e. A topological matrix. A, relating the incidence of element ends 

on the nodes, is now generated as follows. Each of its elements is in 

actuality a 6 x 6 matrix, 0 or I, 0 representing a 6 x 6 null matrix and 

I representing a 6 x 6 identity matrix. The topological matrix will be 

T 

generated in its transposed form, namely A , since both forms are needed 

T 

in subsequent calculations. A will be 6N x 6M in size. The following 
rules obtain: 



ij 



will be I if end j is incident on node i. 



A T . , 
i* j 



will be 0 if end j is not incident on node i. 

I and 0 are 6x6 matrices as defined above. To permit doing this, a 
certain order must be established both for the nodes and the element 
ends, and this order will have to be adhered to for subsequent assemblages 
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and calculations. It will be helpful if the low numbered nodes (i.e., 

T 

the upper rows of the matrix A ) are the boundary nodes as this will 

T 

minimize subsequent reordering. As an example, A is shown for the 
simple system given below. 




Nodes are numbered, element ends are lettered. Nodes 1 and 2 are boundary 
nodes, nodes 3 and 4 are primary nodes. 



abcdefgh 



4 

01 



Row of all element ends 



10000000 
00000001 
01101000 
000I0I10_ 

Column of nodes, note the boundary nodes listed first 
T 

A is the quantity within the brackets. The row of element ends and 
the column of nodes are shown only to indicate how the matrix is generat- 
ed. 

T T 

As implied herein, A is [A ] 



f. All the primary element dynamic stiffness matrices are now 

assembled into a large matrix called the primitive stiffness matrix, K . 

T 

It is necessary that the order of assembly correspond to that of A . 
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2.5 



Defining {fa F b • • •} as Fg (subscript E referring to ends of elements) 



and 



{*a °b *•*} 



as D then 
E 



F_, = K D„ 

E p E 



2.6 



where is the large square matrix in equation 2.5. 



g. The sum of the forces exerted on any node by the element ends in- 
cident upon it must be zero (except for boundary nodes where the constrain- 
ing forces are the "closing" forces) and the displacements of all the ends 

incident upon any given node must be identical. The topological matrices 
T 

A and A permit these conditions to be expressed mathematically. Using 
the subscript N to refer to nodes* just as subscript E refers to ends 
of elements. 



and 



f n ■ A f e 



D_ = A D„. 
E N 



F„ is a column vector of the forces exerted on the nodes and D„ is a 
N N 

column vector of the node displacements. Therefore: 



N 



T 

= A F. 



- A < K P V 

= A 1 K p A D„. 



16 



we then have the re- 



T 

Defining A ^ A as the node stiffness matrix, K^, 
lationship between node forces and displacements: 



N 






2.7 



h. F„ is a column of N 6x1 matrices and all but B of these will 
have all the elements identically equal to zero. F^, and are now 
reordered so that the B uppermost 6x1 matrices of F N and D N are the force 
and displacement vectors respectively of the B boundary nodes. This step 
is unnecessary if the proper order has been previously established as in- 
dicated in 2.e. Partitioning these matrices now yields: 



— 




— 








6B F 
1 NB 


= 


Sll 


^12 




d nb 


6(N-B) q 




^21 


*Sl22 




D N2 


- 










— — 



i. Since boundary constraints are known in the local coordinate 
systems of the various boundary nodes we now transform back to local 
coordinates. Up to this point all reordering has been done using entire 
6x1 or 6x6 submatrices. This has been done in order to preserve the in- 
tegrity of the partitions of the state vector so that we can make the 
transformation back to local coordinates by using the same transforma- 
tion matrices generated earlier to transform the system to global coordi- 
nates. Now with F and D transformed to local coordinates we may find, 

No Ni5 

upon examination of the boundary conditions, that some of the elements of 
Fjjg are also equal to zero. Reordering the matrices again, this time 
element by element instead of by submatrices, so that the non-zero ele- 
ments of F and the corresponding zero elements of D^. occupy the lead- 
ing positions in their respective vectors, and then partitioning the 
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matrices we find: 



F. 



NB 



* 

Sll 


* 

*H12 


* 


* 




^22 









0 




* 




D N2 



2.8 



* 

is a Pxl vector where P is the number of non-zero boundary condi - 
NB 

tions and is less than or equal to 6B. The asterisk is used to in- 
dicate local coordinates and the reordering indicated in the text. 

From equation 2.8 it can be seen that: 







* 








* 




* 


0 




^21 

* 




0 


+ 


*N22 




°N2 



Since the elements of are not zero, and are independent it appears 



that K ^22 must be singular, i.e.. 



**22 



2.9 



* 






*Sl22 


as the frequency determinant, and calling it 





Defining 

for convenience, we can see that it will be of order (6N-P) x (6N-P). 



Since some of the elements of 



M 



are frequency dependent it can be 



seen that there will be some values of frequency for which the frequency 
determinant becomes zero. These frequencies will be the natural fre- 
quencies of the system. We are unable to solve for these frequencies 
directly since it would not be feasible to obtain a literal expression, 
ultimately expressible as a polynomial (or worse) in the unknown fre- 
quencies. Therefore we calculate the value of the frequency determin- 
ant for many different values of frequency (using a high speed digital 
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computer) and plot the value of the determinant versus the frequency. 

From this plot we should be able to pick, at least approximately, the 
natural frequencies. It is then possible to refine these answers by 
making additional calculations over much smaller intervals and plotting 
these results on an expanded scale. If a lumped mass model of the system 
has been used throughout, there will be a finite number of natural fre- 
quencies. If however, a distributed mass model has been used then there 
will be an infinite number of natural frequencies. 

j. It can easily be seen that the sizes of the matrices involved 

in this method may become quite large. The largest single matrix will 

be K , which is 6M x 6M, and even for a system of only fifteen primary 
P 

elements, this matrix would be 180 x 180. A matrix this size has 32,400 
elements and if stored in the core memory of a 32K computer would leave 
only 368 cells for other variables, the monitor, the resident program, 
etc. This is far too few cells for these purposes. If double precision 
arithmetic were to be used then twice as many cells would be required 
and the matrix could not be stored in core memory. 
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CHAPTER III 



REDUCTION OF SUB SYSTEMS 



3.1 Purpose 

It is the purpose of this thesis to develop and exhibit methods by 
which the mathematical manipulations may be simplified or partitioned 
with the particular purpose of keeping well within the core memory capa- 
city of the usual large engineering computer. Matrices too large for the 
core can be handled with the use of auxiliary memories or intermediate 
tapes; but this slows computations by several orders of magnitude and 
is to be avoided if possible, particularly in a problem of this type 
where the solutions are gained by successive iterations. 



3.2 Theory 

The basic propositions which will be demonstrated are; (1) that any 
ganglion or sub-system can be represented by an equivalent stiffness mat- 
rix of size 6B' x 6B' where B' is the number of nodes at which the gan- 
glion is connected to the overall system, and (2) the size of the overall 
system frequency determinant can be reduced so that it is, at most, 12B 
x 12B where B is the number of boundary nodes of the system. 

Start by taking any ganglion of the system that does not include a 
primary element connected to a boundary node, and treat it as if it were 
a complete system having B' boundary nodes, N' nodes, and M'/2 primary 
elements. Going through the manipulations indicated in part II the 
equation representing this ganglion eventually may be reduced to the form: 



V " V V 



(cf. equation 2.7) 
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Then by ordering the matrices so that the B' boundary node vectors appear 
first in F N and D N and partitioning we obtain: 



f nb' 


= 


^11 


^N12 




°NB' 


_ F N2_ 




^21 


^22 




_ D N2_ 



3.1.1 



Now it must be understood that there is a difference between this case 
and the solution of the overall problem, the difference being that while 
the elements of F^ are identically equal to zero, the elements of 
D , are not. In other words, for the overall system the existence of 
a boundary force indicates a rigid constraint component;? (i.e., no dis- 
placement) whereas for the sub-system, there is no such restriction. 
Proceeding, 



f nb' 




*H11 °NB' + 


CM 

d 3 

CM 

r— 1 

uP 


0 


S3 


\n d nb' + 


CM 

25 

Q 

CM 

CM 

J 3 


Solving the 


second 


of these for D„„ we 
Nz 


get 


°N2 " -K N22 ^N21 °NB' 

Substituting it into the first equation yields 


F . » 
NB' 


[>1 


^N12 ^22 ^ 21 ] 


d nb' 



k n- 



NB' 



3.1.2 



where F^, and D^g, are 6B' x l and K^, is 6B' x 6B' . K^, is defined by 

this relation. 

What has been generated might be thought of as a B' -ended pseudo- 
element to replace the ganglion and retain its identical characteris- 
tics. The pseudo-element is represented by a stiffness matrix, but note 
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that this is no longer the 12 x 12 stiffness matrix of a simple element 
but a 6 b* x 6B' matrix. 

This procedure may be repeated for other ganglia which are separate 
from those considered previously, or a pseudo-element and its associated 
stiffness matrix (or indeed several such elements) may be used as a build 
ing block in constructing a larger ganglion. Each pseudo-element calcula 
tion eliminates N'-B' nodes from the sub-system. This leads to the con- 
clusion that the most "efficient" ganglion to choose is the one with the 
fewest boundary nodes for the most internal nodes. 

With the various ganglia reduced to pseudo-elements, it is then 
possible to proceed as in Chapter II. A primitive stiffness matrix is 
assembled and pre-and post-multiplied by the topological matrices to 
yield a node stiffness matrix. This is reordered, transformed, and re- 
ordered again, the frequency determinant is isolated, and the natural 
frequencies are found by an iterative technique as described in section 
2.i. What actually has been done is to build up to a solution step by 
step, eliminating as many nodes as possible along the way. It is obvious 
that for the final step one can have the B boundary nodes of the system 
each connected by a primary element to one of the B' boundary nodes (B'^ 
B) of a pseudo-element representing the entire remainder of the system, 
i.e., the entire system less the elements connecting points B with points 
B' . 

It can b« shown that this is so by considering the system as a 
whole and then reducing it, rather than by starting with a small sub- 
system and building up to the complete system. The first ganglion can 
then be chosen to be everything except those primary elements that have 
one end incident on a boundary node. If such a ganglion can not be 
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chosen, then there is more than one structurally independent system in- 
volved and they must be solved separately. This statement is not as 
trivial as it may seem to be at first reading. It is entirely possible 
that a piping network that would be considered as one system by reason 
of its application or fluid flow could actually be several independent 
structural systems. This could be caused by the pipe being rigidly 
anchored at one or more points thus effectively isolating the pipe on 
one side of the anchor from the effects of any change on the other side. 

Although many of the boundary nodes of a piping system may be topolo- 
gically identical (for example all nodes at which an element is rigidly 
connected to "ground") it is necessary to consider each as a separate 
entity because of the necessity of using different local coordinate 
systems and then different transformation matrices. Thus if two or 
more elements terminate on what is geometrically the same boundary node, 
they must be considered as terminating on topologically separate nodes. 
Therefore, there will be only one primary element incident on each bound- 
ary node. (See definition, page 7 ). Keeping this in mind it can be 

seen that the number of nodes connected to boundary nodes by primary 
elements must be less than, or at the most equal to B. Thus for the 
reduced system N is less than, or at most equal to 2B. Taking the 
limiting case, N = 2B, the node stiffness matrix will be 12B x 12B 

and the frequency determinant will be (12B-P) x (12B-P) where P lies in 
the range 0 ^ P^6 b. If all the boundary nodes represent complete con- 
straints in all degrees of freedom P will be 6B, the other limit being 
no constraint on the boundaries and P equaling zero. Thus the order of 
the frequency determinant will be between 6B x 6B and 12B x 12B. The 



23 



restriction of having only one primary element incident on a boundary 
node would not be necessary for the solution of a mathematically similar 
problem where physical location and metric geometry are of no importance, 
for example, an electrical network. 

3.3 Reduction of Parallel Elements 

To provide a tool for future computational work, and to demonstrate 
the use of topological matrices it will be shown that parallel elements, 
i.e., primary elements connecting the same pair of primary nodes, can be 
reduced to an equivalent element. The stiffness matrix of this equi- 
valent element will be the sum of the stiffness matrices of the various 
component primary elements. This is not a demonstration of the reduc- 
tion method of section 3.2 as there are no nodes eliminated. The re- 
duction of two parallel elements into one equivalent element is shown 
here. The logical extension is then made that this equivalent element 
can be taken parallel with another primary element and these may be re- 
duced to a second equivalent element. As many parallel elements as 
required can be handled in this manner. 

Take the case of two parallel elements. Fig. 3.3.1. 

a ^ \b 

1 c J> 2 Fig. 3.3.1 

cX - — ''a 

1 and 2 are nodes, and a,b,c, and d are element ends. 



" U ba 
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from which 


F K 
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ab 


D a,b 


“ U dc 
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F ^ 
c,d 


a 

o 

CL 


D c,d 



After transforming these to some global coordinate system the primitive 
stiffness matrix can be assembled, viz.: 
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3 . 3.2 
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it has been shown that 




3 . 3.3 



Thus, the equivalent stiffness matrix for two parallel elements is the 
sum of the two individual stiffness matrices. From the above we can see 
that the stiffness matrix for n parallel elements would be the sum of the 
n individual stiffness matrices. 

This intuitively attractive result could probably have been shown 
without resorting to the topological matrices. The preceding deriva- 
tion, however, puts into a definite mathematical form what would have to 
be understood in another derivation. That is, the sum of the forces on 
the node equals the sum of the forces at the element ends and the dis- 
placements of all element ends at any node are identically equal. 

An equivalent element with its own stiffness matrix has now been 
generated, and may be used as one of the building blocks of the system. 
However, in a manner equivalent to reversing what we did to get a stiff- 
ness matrix from a transfer matrix in section 2.1 we can get, if desired, 
a transfer matrix for this equivalent element, in whichever local co- 
ordinate system is convenient. If the system is so configured that after 
the reduction of the parallel elements the equivalent element is then one 
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of two or more elements of a primary path, the transfer matrix may then 
be used in the generation of a primary element including the equivalent 
element. 
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CHAPTER IV 



COMPUTER PROGRAM PLANNING 



4.1 Discussion 

It 'Should be obvious by now that the solution of any but the simplest 
problems by these methods will require the use of a high speed digital 
computer, and that the use of intermediate tapes or auxiliary memory will 
most likely be required. 

The primary object of all the calculating is to find some or all of 
the real values of frequency which will make the frequency determinant go 
to zero. The corresponding "mode shapes" can easily be computed by a 
repeated pass using the correct frequency. Since the computer can only 
calculate with numbers, it will not be possible to set up an equation for 
the determinant and solve it exactly. It will not even be possible to set 
up an equation and try to solve it by substituting very many values of 
frequency. Rather, all calculations must be done numerically using assum- 
ed values of frequency. The frequencies which make the frequency deter- 
minant equal to zero must be found by an iterative process. This means 
that virtually all the calculations must be redone for each assumed fre- 
quency. 

The writing of a computer program to do the work involved is not a 
simple matter. The ideal program should be written in a widely used en- 
gineering language (such as FORTRAN) and arranged to allow data to be 
presented in a most expeditious manner. A good example of this type of 
program is "STRESS" by Fenves et. al. [4] [5], "STRESS* 1 solves static 
force-displacement problems by use of a stiffness matrix method. 
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No attempt has been made to write a program for this thesis. How- 
ever, the outline of such a program is shown in Fig. 4.1.1. A method for 
systematically approaching this problem is given in section 4.2. 

The reader interested in actually programming this type of problem 
is advised of the existance of a FORTRAN program written by Fink [6] in 
1964. Program VIPIPE is severely restricted in the size and topological 
configuration of the systems it can handle. It was written to determine 
in-plane and out-of-plane vibrations of a two dimensional piping system 
and does not make use of topological matrices. It would be a good start- 
ing point for the generation of a more sophisticated program and many of 
its subroutines, in particular those for calculating transfer matrices, 
could probably be put to good use. 
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Figure 4 . 1.1 
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4. 2 Expansion of Computational Concepts 



The purpose of this section is to develop to a greater extent the 
approach to a computer solution for a general problem of this category. 
The basic flow diagram for a computer program shown in Fig. 4.1.1 will 
be used as a rough guide. 

The following assumptions are made: 

a. A computer memory of about 30K available cells is to be 
utilized. 

b. Intermediate tapes rather than drums or discs will be 
used. This assumption puts some constraints on the 
intermediate storage of information as binary tapes cannot 
be randomly accessed whereas drums or discs can. 

c. Most of the calculating will be done with double precision 
arithmetic. 

d. The system to be solved may be of any size. 

The following systematic approach is recommended. 

1. ) Display the system topologically showing only the boundary 
and primary nodes and the primary elements connecting them. Number the 
nodes, using the lowest numbers for the boundary nodes, and number the 
primary elements. 

2. ) Determine the way in which the system will be divided into 
ganglia and subsequently recombined. A procedure for making cots for 

an "optimum" reduction is beyond the scope of this work, and indeed per- 
haps the idea of an optimal reduction is frivolous. However, certain 
operations should be readily apparent, and the first of these is the 
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combining of parallel elements, where present. From this point forward 
it is an empirical process. One must use his judgement in working to- 
ward replacing the entire system, except for the boundary nodes and the 
primary elements incident on the boundary nodes, with a single pseudo- 
element. For reasons that will be discussed later, it is best to set 
up two or more different sets of ganglia so that the problem may be solved 
more than once and the answers compared. 

3. ) Display the system so that the geometry will be readily avail- 
able for computer programming. An isometric line drawing, (not neces- 
sarily to scale, but at least reasonably proportional to the actual layout) 
will probably be adequate. 

4. ) Analyze the boundary node constraints. It is necessary that 
each node be either unrestrained or rigidly restrained in each of its 
six degrees of freedom. If this is not the case then an imaginary ele- 
ment with characteristics equivalent to the desired boundary constraints 
must be added between the element end and the boundary node to satisfy 
the above condition. 

5. ) Analyze each other node. The assumption made in this method 
is that elements are rigidly connected together at the nodes. If this 
is not the case, an imaginary spring element of desired characteristics 
must be included wherever required.* 

6. ) Establish an origin of coordinates and an orthogonal set of 
global coordinates. 

*This is related to the difficult problem of releases. See reference 
[7]. 
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7. ) Determine and tabulate for each element the following data: 

a. Type of element 

1. straight pipe 

2. arc of circle pipe 

3. "imaginary" element 

4. non-pipe (hanger) 

b. weight per unit length of pipe 

c. weight per unit length of fluid and insulation 

d. pipe diameter 

e. wall thickness 

f. E (modulus of elasticity) 

g. G (modulus of elasticity in shear) 

h. coordinates of ends and of center of curvature (if an 
arc) 

i. radius of curvature 

8. ) Read in the data for each element; make whatever calculations 
are necessary to convert the raw data into an efficient form for future 
calculating; and store all this information on tape. All this data is 
independent of frequency, and must be used in each iteration to calcu- 
late transfer matrices. It is thus vital that the information be stored 
in the order that it is to be used so that the tape may be used in the 
most expeditious manner. Let us assume that we have a set of data cards 
carrying all the information for all the elements of each primary ele- 
ment. We also have a larger set composed of all these sets, so that the 
data is on cards for all elements of the system. Now go back to 2) and 
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pick the first ganglion to be reduced. Decide on the order in which the 
primitive stiffness matrix is to be assembled. (This will also be a 
good time to determine the way in which the node stiffness matrix is to 
be ordered and to generate the topological matrices for this ganglion, 
though these items will not be needed for a while.) Assemble the data 
cards for each primary element in order from one end to the other. 

Remember that the coordinates of the first end of the first element will 
determine the local coordinate system in which that particular primary 
element will be calculated. For each primary element set there should 
be an initial card indicating the numbering of elements in the primary 
element. Assemble all the primary element card sets in the order in 
which the primitive stiffness matrix is to be assembled. 

9. ) Choose the next reduction to be made. If it is independent of 
the primary elements already assembled for reading, proceed as above. If 
it contains the ganglion (or several ganglia) assembled for reading, 
proceed as above only for the primary elements not already used. If it 
is composed entirely of elements whose cards are already assembled, no 
action is required at this time. ( It would still be good at this time 
to set up the order for the node stiffness matrix and generate the topo- 
logical matrices for this ganglion) . 

10. ) When the ordering is complete the element data is ready to be 
read in. A flag indicating that no more element data is to be read should 
be added. The number of elements in each primary element are stored in 

an array to facilitate retrieval of information from the tape. The decis- 
ion on exactly how to handle data, particularly in reference to the amount 
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of data to be read in and transferred to tape at one time must be left to 
the programmer. It is faster to accumulate data in core and transfer it 
all to tape with one WRITE TAPE statement than to go to tape with each 
individual set as it is calculated. It should be noted that much economy 
in data card preparation can be obtained by eliminating the necessity for 
duplicating information that is identical for many elements, i.e.. Young's 
modulus, diameter, etc. 

11. ) Now starting with the first ganglion to be reduced, read in 

T 

the appropriate topological matrix (A or A ) and store it. It will not 

T 

be necessary to have both A and A in memory since by programming we can 
utilize the stored matrix either as itself or as its transpose. Since all 

or j>0 great storage 
economy can be obtained by representing these 6x6 elements by a Boolean 
1 or 0 respectively and storing thirty-two of these per core memory word.* 
It will probably not be necessary to go out on tape with these matrices. 

At this point all the data that will be constant from iteration to 
iteration should be either on tape or in core. 

12. ) With the assumption of an initial frequency the first iteration 
may now be started. From this point on, all calculations will be done 
with double precision arithmetic. Start by taking the first ganglion as 
previously ordered, and read back into the core all the data relevant to 
the elements of this ganglion. By the use of sub-routines the stiffness 
matrix for each primary element of the ganglion is calculated. It will 

*The number of logical elements that may be stored in one word varies 
from computer to computer. The number thirty-two used here applies to 
the CDC 1604. 



elements of the topological matrices are either gl 
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not be necessary to store the numerous transfer, rotation and transforma- 
tion matrices that are generated for if they are generated in the correct 
order and multiplied together then all that need be stored is the pro- 
duct. The stiffness matrix for each primary element will occupy 2 x 12 
x 12 - 288 core storage locations where the 2 results from double pre- 
cision arithmetic. Now although the primitive stiffness matrix, K , is 

P 

shown as a large square matrix, it will be assembled in the computer as 
a three dimensional array. The reason for this is that K p consists only 
of diagonal elements, each of which is of order 12 x 12, and all the other 
elements are zero. By using a three dimensional array we can store only 
the diagonal elements and not waste space storing zeros. The storage re- 
quired for Kp will be 288 x M'/2, and if the stiffness matrices have 

been cleverly stored when they were generated already exists. The 

T 

node stiffness matrix is now generated from a Ak^A. Note that the 
column vectors of node forces and displacements are never actually used 
in the solution. From we generate K^, as shown in equations 3.1.1 
and 3.1.2. The order of this matrix will be 6B' x 6 b'; and the number of 
memory words necessary to hold it will be 2 x 6B' x 6B'. As before, the 
2 represents the requirement of two words per item in double precision 
arithmetic. 

Assuming that this reduced ganglion will be used as part of the 
next ganglion to be reduced, we will plan to leave K^, in core memory. 

If it is not to be used for a while, and if space is at a premium for 
other calculations, it may be better to store K^, on tape. 

13.) By continuing in a similar manner eventually we will arrive 
at the largest ganglion, namely, that one including everything except 
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the primary elements incident upon the boundary nodes. The order of the 
resultant stiffness matrix for this ganglion will be 6B x 6B. Using 
this stiffness matrix, and the stiffness matrices of the remaining 
primary elements we now set up the primitive stiffness matrix for the 
entire system and by pre- and post-multiplying it by the topological 
matrices obtain the node stiffness matrix for the entire system. Since 
we have assembled these matrices carefully, the reordering described in 
section 2.h will not be necessary. The transformation back to local co- 
ordinates is made for the boundary nodes. With the boundary nodes re- 
presented in local coordinates it is now possible to make the second re- 
ordering as described in section 2.i. The frequency determinant is now 
isolated and its value calculated. This value and the corresponding 
value of frequency at which the calculations have been made are saved. 

14.) A new value of frequency is now assumed and the above process 
repeated. Subsequent values of the determinant may be used in conjunc- 
tion with the values found previously to predict the next frequency to be 
tried in order that the calculations may converge on the desired natural 
frequency. Alternatively, calculations may be made at constant intervals 
and the results plotted. 

4.3 Computational Difficulties 

There is one predominant problem that keeps appearing, if only 
implicitly, in work with matrices. Any time matrix equations are mani- 
pulated in such a way that the inverse of some matrix is required, the 
assumption must be made that the matrix to be inverted is non-singular. 
Matrices are inverted many times in the course of a solution by the methods 
discussed herein, and the possibilities of having singular matrices occur 
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occasionally are significant and somewhat disturbing. 

In the determination of the stiffness matrix of an element from its 
transfer matrix it is necessary to invert a 6 x 6 matrix as described in 
section 2„b. 

non 

m 
F 

l_ raJ 

from which 

K 

Note that i® the 6x6 matrix to be inverted. As some or all of the 
elements of are frequency dependent (all other parameters are depen- 
dent upon physical constants and are thus fixed) it can be seen that there 
will be some values of frequency for which the determinant of Uj^ will 
vanish. These frequencies correspond to the natural frequencies of the 
element with both ends "built-in" or cantilevered. In matrix form, from 
2 . 2 , 
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12 



= 0. Only the real values of frequency that satisfy 
this relationship are of interest. 

What this means in terms of programming is that whenever the assum- 
ed value of frequency at which a series of calculations is being made is 
the same as (or close to) one of the natural frequencies of one of the 
primary elements it will not be possible to generate the stiffness matrix 
for that element. It will be necessary to assume a new frequency near the 
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unacceptable one and to try to calculate to a conclusion again. 

A similar problem may arise in the calculation of the equivalent 
stiffness matrix for a ganglion. Assuming that stiffness matrices have 
been satisfactory calculated for each primary element of the ganglion 
as shown in the arithmetical manipulations of equations 3.1.1, it is 
still necessary to invert a matrix and this matrix, also being made up of 
elements dependent upon the frequency, may become singular. The fre- 
quencies at which this matrix becomes singular correspond to the natural 
frequencies of a system equivalent to the ganglion but with the boundary 
nodes rigidly fixed. 

The inability to calculate to a conclusion due to an inability to 
invert part of the node stiffness matrix of a ganglion can be avoided, 
in most cases, by breaking the system up differently and thus using dif- 
ferent ganglia. It is for this reason that in section 4.2 it is recommend 
ed that any problem of this type be solved using several different combina 
t ions of ganglia. Any reasonably sophisticated computer program will, of 
course, provide some sort of diagnostic message and proceed to calculate 
at a different value of frequency whenever some matrix cannot be inverted. 

The problem of a singular matrix affecting the calculations should 
not be too serious at the low frequency end of the vibrational spectrum. 

It can reasonably be expected that the natural frequencies of a ganglion 
will be considerably greater than the lowest natural frequency of the 
overall system. It is not possible to forecast at what frequency pro- 
blems will actually begin to occur. What is to be done is just to pro- 
ceed with the calculations, avoiding those frequencies at which calcula- 
tions cannot be made, but calculating as closely as possible on both 
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sides of such frequencies. 

There are other problems involved in the creation of a computer pro- 
gram capable of handling large linear systems by matrix methods. There 
is, for example, the possibility that even after all reductions have been 
made the frequency determinant will still be too large for core memory. 
There will also be the need for making a number of tests during the course 
of a program run, both for debugging purposes and to provide a hedge 
against gross errors. (An example of this would be the testing of the 
value of mass of an element to insure that it is positive). It is re- 
commended therefore that the ’’team" for developing this program consist 
of at least a mechanical engineer well versed in vibrations and a com- 
puter programmer with experience in handling large systems and matrices. 
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